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Abstract 

A numerical hydrodynamical model for the evolution of spherically symmetric 
collapsing clouds, designed for the calculation of the thermal structure of these ob¬ 
jects in both the prestellar and protostellar stages of their evolution, is presented. 
Distinctive features of the model include the possibility of independently describing 
the temperatures of the gas and dust, which is extremely important when calcu¬ 
lating the thermal structure of prestellar and protostellar clouds, and the account 
of the radiation flux from the central protostar. This model is used to compare 
the theoretical density and temperature distributions with observations for nearby 
sites of star formation obtained with the Herschel Space Observatory. Application 
of the diffusion approximation with a flux limiter describes well the radial density 
and temperature distributions in protostellar clouds. However, significant differ¬ 
ences between the model and observational density profiles were found for prestellar 
stages, suggesting the presence of appreciable deviations from equilibrium in the 
prestellar clouds. An approximate method for calculating the thermal structure of 
a cloud based on the adaptive r-approximation is presented. Application of the 
r-approximation yields good agreement with the diffusion approximation for the 
prestellar phase, but produces appreciable discrepancies for the protostellar phase, 
when the thermal structure of the accreting envelope is determined by the radiation 
of the protostar. 

1 Introduction 

A model for calculating the thermal structure of a collapsing protostellar cloud was 
presented in [1, 2]. This model is fast enough for use in hydrodynamic calculations, 
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and is also sufficiently accurate to use the modeling results to calculate theoretical 
spectra that can be compared with observations. The key feature of the model is 
separate treatment of dust and gas temperatures, since these two temperatures may 
differ significantly in early phases of the evolution of the clouds, and also in their 
outer regions. Different heating and cooling mechanisms are considered for the gas 
and dust, while they exchange energy through collisions. Since the main source of 
heating (cooling) of the dust is the absorption (emission) of radiation, the method 
used to calculate the radiation transfer is an important component of the model. 

The basic idea of the radiative transfer computation method is to divide the 
overall frequency range into low-frequency (infrared) and high frequency (ultra¬ 
violet) parts. The high-frequency part is diluted by interstellar radiation. The 
intrinsic radiation of the dust can be neglected in this part of the spectrum, and 
the absorption (and scattering) of external radiation alone can be considered. In 
the low-frequency part of the spectrum, it is necessary to consider both absorption 
and thermal emission of the ambient medium. This division of the spectral range 
into two parts enables the use of different appropriate approximate methods in each. 
Modeling the UV part of the spectrum reduces to calculating the mean intensity 
of the UV radiation via direct integration of the radiative transfer equation along 
selected directions. The calculation of the intensity of the IR radiation is based on 
solving a system of moment equations that represent the diffusion approximation. 
Thus, the model contains four interacting components: gas, dust, and the IR and 
UV radiation. 

An axially symmetric realization of this thermal model was used to calculate 
the early evolution of a contracting, magnetized protostellar cloud [1]. Modeling 
of later stages, up to the formation of the first hydrostatic core, was performed in 
a spherically symmetric approximation in [2]. In our present study, we describe 
a modification of the model designed for the calculation of the accretion phase in 
the evolution of a protostellar cloud following the formation of the protostar. We 
present the results of our calculations of the structure of a spherically symmetric, 
contracting protostellar cloud and accreting envelope. To simplify the computation 
of the thermal structure, we propose the adaptive r-approximation method and 
compare it with the original method. We also compare the model and observed 
distributions of the density and temperature in prestellar and protostellar cores of 
molecular clouds. 


2 MODIFICATION OF THE MODEL 

Pavlyuchenkov and Zhilkin [2] investigated the thermal evolution of a protostel¬ 
lar cloud up to the formation of the first hydrostatic core. Evolution of this core 
subsequently leads to the evaporation of dust, dissociation of hydrogen, and, con¬ 
sequently, to the collapse of the first core and the formation of a second hydrostatic 
core, which is essentially the protostar [3-5]. Simulation of these stages requires 
that additional physical processes be taken into account in the model for the proto¬ 
stellar cloud. The timescales for these stages are extremely short compared to the 
dynamical timescale of the accreting envelope. This makes hydrostatic modeling 
of the core and accreting envelope using a single model extremely difficult. More¬ 
over, multi-dimensional models must be used in studies of hydrostatic cores, since 
their formation and evolution are closely associated with the formation of accretion 
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(protoplanetary) disks and outflows from them. However, in addition to modeling 
the evolution of a protostar and its surroundings (the protostellar disk), studies of 
the structure of the extended accreting envelope (including its thermal structure) 
are of interest. Accretion from the envelope supplies matter to the protostar and 
protostellar disk, determining their evolution. In turn, the observational identifi¬ 
cation of protostellar accreting cores depends on assumptions about their thermal 
structure, which requires appropriate models. 

Historically, the simplest models for prestellar and protostellar envelopes were 
based on semianalytical self-similar solutions of the hydrodynamical equations in 
isothermal or polytropic approximations [7, 8], with possible modifications to take 
into account rotation in the cloud [9], as well as numerical hydrodynamical model¬ 
ing taking into account the finite size of the cloud [10-12], Later, numerical models 
were used to consider the nonisothermal structure of the envelopes more accu¬ 
rately, via more accurate solution of the transfer equation taking into account the 
frequency dependence of the absorption coefficient [13], and applying the diffusion 
approximation with a flux limiter and either mean [14] or frequency-dependent [15] 
opacities. Modeling of protostars and protostellar envelopes has been carried out 
by a number of research groups (see, e.g., [6]). However, most models still assume 
that the temperatures of the gas and dust are equal, since frequent collisions should 
establish thermal equilibrium between these components. As is noted above, how¬ 
ever, this assumption may not be valid at some stages of the evolution. In some 
models, the dust and gas were treated separately by postprocessing the gas density 
distribution in the envelope using Monte-Carlo methods for the radiative trans¬ 
fer [16]. Finally, most studies have focused on early stages in the evolution of the 
first hydrostatic core and protostar [17]. Here, we consider a collapsing cloud and 
describe a modification of our model [2] designed for the calculation of the thermal 
structure of an accreting envelope. 

The model of [2] is based on a numerical solution of the hydrodynamical and 
radiative transfer equations. The continuity equation, equation of motion, and the 
equation describing the variations of the thermal energy of the gas E g have the 
form 


d_E g 

dt 


dp 

dt 


+ V • (pv) = 0, 


<9v . Vp 

_ + (v . V)v = __ + g . 


+ V • (E s v) + pV • v = T g — A g — A gd . 


(1) 

( 2 ) 

(3) 


When the velocity of the medium is low, the equations for the thermal energy 
of the dust Ed and the energy of the IR radiation E r have the form 


— — <7p(aT d — E t ) + Agd + S, (4) 

^ = v- (^ v ^ r ) +^p(ar d 4 -^)- (5) 

In these equations, p is the density, v the velocity, p the pressure, g the gravi¬ 
tational acceleration, T g the gas heating function, A g the gas cooling function, A gd 
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the rate of exchange of thermal energy between the gas and dust, S the rate of dust 
heating by UV radiation, the dust temperature, ap and ap coefficients that 
depend on the opacity of the medium, and a the radiative constant. A detailed 
description of the system of equation and the solution method is given in [2]. Here, 
we will discuss our modification of this model. 

The modification was carried out in the framework of the well-known “sink cell” 
formalism (see, e.g., [18]). In this approach, the central region of the cloud (the 
sink cell) is excluded from the simulation when a high density is achieved there. 
This makes it possible to exclude from the modeling the rapid evolution of this 
region with numerous physical processes leading to the formation of the central 
star. However, the influence of this region on the surrounding cloud is taken into 
account, as well as the influence of the cloud on the inner region. This interrelation 
has both a dynamical and thermal nature. Thus, the modification of the model 
affects the hydrodynamical method and the method used to compute the thermal 
structure. The modification of the hydrodynamical method depends on the type 
of method used (Eulerian, Lagrangian, smoothed particle hydrodynamics method, 
etc.). 

We will now describe the changes in the one-dimensional Lagrangian method 
used. Let r ce ii be the radius of the inner sink cell. As the cloud contracts, the cells 
of the computational domain become denser and move toward the center. We will 
suppose that a cell is excluded from the calculations (it is absorbed by the sink 
cell) as soon as the right-hand border of the cell crosses the radius r ce n. The next 
cell then becomes the boundary cell, and the cells are renumbered such that the 
new boundary cell is assigned the subscript “1”. Thus, the left-hand boundary of 
the new boundary cell r\ will be located to the left of r ce u, while its right-hand 
boundary r 2 is located to the right of r ce u. We assume that the pressure gradient 
vanishes at the left border (ri): 


dp 

dr 


= 0 . 


n 


( 6 ) 


This condition means that the matter flows freely (by inertia) into the inner 
zone under the action of gravity of the sink cell alone. 

The modification of the method used to compute the thermal structure involves 
two operations: taking into account the sink cell and modifying the diffusion op¬ 
erator. Taking into account the sink cell reduces to formulating a new boundary 
condition at the inner boundary of the computational domain n, where the IR flux 
from the sink cell L ce \\ is determined by the sum of the photospheric and accretion 
luminosities of the young star: 


-7'cell — 1-star T Lacc- (7) 

We calculated the photospheric luminosity using the stellar model of [19, 20], 
where the luminosity L star and the radius of the young star R star are computed as 
functions of the instantaneous M star and age 4tar- Assuming equipartition between 
the thermal and kinetic (mechanical) energy of the accreting mattei0, the accretion 
luminosity is given by 

^The mechanical energy does not contribute to the radiation and is mostly carried away together with 
polar jets and outflows of matter. 
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lGAfM star 


(8) 


-rt'Star 


where 


M = AirrfpiVi 


(9) 


is the accretion rate through the inner boundary of the computational domain. 
Here, p\ is the density in the first (boundary) cell, and v\ is the radial velocity at 
the left-hand border of the first cell. 

When calculating the accretion luminosity using these formulas, we assumed 
that all the matter crossing the boundary of the computational domain is in¬ 
stantaneously accreted onto the star, and the all thermal energy of the accreted 
matter is converted into IR radiation; this corresponds to the so-called cold ac¬ 
cretion limit, and is valid if the accretion rate onto the protostar is not too high, 
M < HR 5 Af 0 /yr [21], Since the coordinate in the Lagrangian method is the 
mass inside the corresponding radius, the mass inside the absorbing zone can be 
calculated directly from the coordinates of the boundary cell: 


( 10 ) 


Afstar — 47TQ1, 


where q\ is mass (Lagrangian) coordinate. 

Finite-difference approximation of the diffusion operator in the spherically- 
symmetric case 




( 11 ) 


A E 


where 



( 12 ) 

(13) 

(14) 

(15) 


Ri ~ r i-l/2 r i+l/2i 
An = r i+ 1 /2 - ri_ 1/2 , 
Ar i+ 1/2 = r i+ 1 - n. 


Taking into account the new boundary condition, the finite-difference approxi¬ 
mation of this operator for the inner boundary of the cell takes the form 



(16) 


We introduce the sink cell when the central density becomes n(H 2 ) = 10 14 cm 


The central temperature at this time is close to 1000 K. The radius of the sink cell 
2 Note that there is an error in Eq. (38) in [2]. 
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we used is r ce ii = 50 AU. This is appreciably larger than the radius of the first 
hydrostatic core, but appreciably smaller than the extent of the accreting envelope. 
This choice ensures that the sink cell is larger than the first hydrostatic core. Oth¬ 
erwise, the formation of a quasi-hydrostatic structure at the inner boundary of the 
computational mesh would hinder the accretion of matter onto the sink cell. When 
the sink cell is introduced, all the grid cells inside this zone are eliminated from 
further consideration. The time when the sink cell is introduced is also taken to 
be the time when the star forms, when we start to compute the photospheric and 
accretion luminosity. 

Other numerical criteria for the formation of the protostar are also possible. In 
particular, it was found in [4] that the maximum mass of the first hydrostatic core 
(i.e., the mass at the beginning of the dissociation of molecular hydrogen and the 
subsequent formation of the protostar) was close to 0.05 M q , independent of the 
initial mass of the cloud. This can be used to define the moment of the protostar’s 
formation. However, this value was obtained for nonrotating, spherically symmetric 
clouds. Taking into account the angular momentum of the cloud could lead to a 
reassessment of the maximum mass of the first hydrostatic core [22], 

The Eddington approximation was applied in the equation for the energy of the 
IR radiation in [2]. This approximation assumes that the radiation field is close to 
isotropic, making it possible to express the radiation pressure tensor P r in terms of 

the radiation energy E r in the form P v = —IE r , where / is the unit tensor. The 
diffusion coefficient in (5) is then 


De = 1 /ctr 


c 

ZpdKR 


(17) 


where pd is the dust density, kr the Rosseland mean opacity, and c the speed of 
light. However, the Eddington approximation yields large errors when the radiation 
propagates from a localized source in an optically thin medium (in the so-called 
streaming regime). The accretion stage in the evolution of a protostellar cloud 
proceeds in precisely this regime, with the protostar playing the role of the localized 
source and the accreting envelope being optically thin in the IR. Therefore, in place 
of the diffusion coefficient in the Eddington approximation, we used a flux limiter 
(so-called Flux Limited Diffusion, FLD) [23]. 

If the frequency-averaged scattering coefficient is substantially lower than the 
absorption coefficient (as is true in the IR), the diffusion coefficient is given by 


where 


Dfld 


c aT d 
PdKp E r 


2 + / Vgr 

6 + 3 f + p' fMKpaT^ 


(18) 


(19) 


where up is the Planck mean absorption coefficient. In the iterative scheme 
used to derive T^ +1 and E" +1 , the values of -Dfld were computed using the dust 
temperature TT d and the energy of the IR radiation E™ from the previous time 
step. The FLD approximation operates well in both the streaming regime and the 
optically thick regime, when Dfld ~ Dp. Thus, the evolution of the accreting 
envelope and the preceding stage of the formation of the hydrostatic core should 
both be computed correctly in the FLD approximation. 
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Figure 1: Distributions of the hydrogen number density (left) and the velocity (right) in 
the protostellar cloud at various times during its evolution. The dash.dot lines show the 
functions n oc r~ 3 / 2 and v oc r -1 / 2 , whose indices are dehned by the self-similar solution 
in the case of spherically symmetric accretion. 

3 EVOLUTION OF THE CLOUD WITH TRAN¬ 
SITION TO THE ACCRETION PHASE 

Figures f and 2 present the results of our modeling of the evolution of a protostellar 
cloud with the initial conditions taken from the model described in [2]. The distri¬ 
butions of the hydrogen number density n(H 2 ), velocity v, and the temperatures 
of the gas (T g ), dust (Td, and IR radiation (T r ) are shown for four times. Let us 
describe the main features of these distributions for each of these times. 

1 . t = 0.0 ffj. The initial protostellar cloud with parameters identical to those 
described in [2]. The distributions of the density and temperature were obtained 
for a cloud with central density n(H 2 ) = 10' 5 cm -3 in hydrostatic and thermal 
equilibrium. The cloud was irradiated by external blackbody radiation with the 
temperature fO 4 K and dilution 10 -14 . To initiate the contraction, the density in 
the hydrostatic cloud was increased by a factor of two. The temperature of the 
dust and gas in the initial configuration are different throughout volume, due to 
the relatively low density of the medium. The thermal structure is governed totally 
by the external radiation. 

2. t = 1.6 tfj. The first hydrostatic core. As a result of the cloud compression, 
the central hydrogen number density has reached 10 14 cm -3 , and the tempera¬ 
ture 1000 K. While the temperatures of the gas and dust in the inner regions 


7 








io ' 1 io u io' icr io° icr io° 

R, AU 


io ' 1 io u io 1 icr io° icr io° 

R, AU 



R, AU R, AU 

Figure 2: Distributions of the gas (T g ), dust (T d ), and radiation (T r ) temperatures at 
various times during the evolution of a protostcllar cloud. The dash-dot line corresponds 
to the temperature distribution of the stellar radiation in the absence of absorption. 









(r < 300 AU) are equal, they remain different in the outer parts of the cloud. 
The temperature of the IR radiation is close to the temperature of the gas and 
dust inside r < 300 AU, indicating that this region is opaque to its own thermal 
radiation. The quasi-hydrostatic equilibrium of the inner region can be inferred 
from the velocity profile: the contraction rate is close to zero at r < 10 AU. This 
epoch corresponds to the time when the sink cell was introduced. 

3. t = 1.7 iff. The time 0.1 iff after the formation of the protostar and introduc¬ 
tion of the sink cell. The temperatures of the gas, dust, and radiation are equal 
within r < 1000 AU, and reach 100 K at the inner boundary. At the same time, 
the temperatures in the envelope are close to those in the prestellar phase of the 
evolution. The hydrogen number density and gas velocity at the inner boundary 
are 10 8 cm -3 and 2 km/s. 

4 . t = 2.7 iff. The protostellar phase, corresponding to a developed accretion 
regime. At this time, the mass of the star is 3 M & , its age is 90000 yr, its radius 
11 Rq , the accretion rate 3.6 x 10” 5 Mq/ yr, the photospheric luminosity 64 L@, 
and the accretion luminosity 155 Lq. The distributions of the density and velocity 
are monotonic and are described well by the power laws p oc r” 3 / 2 and v oc r” 1 / 2 , 
whose indices are characteristic for the self-similar solutions in the case of spherically 
symmetric accretion [7] (they are shown by the dash.dot line in Fig. 1). At the 
border of the sink cell, the hydrogen number density is close to 10 8 cm” 3 , the 
accretion velocity has increased to ~ 10 km/s. The temperature of the gas and 
dust in the inner regions is governed by the IR radiation of the protostar, while it 
continues to be determined by the interstellar radiation at the outer boundary of 
the cloud. The dash-dot line in Fig. 2 shows the distribution of the temperature of 
the stellar radiation in the absence of absorption, calculated using the expression 

T(r) = ( cc I r” 1 / 2 . The model distribution is located slightly above this 
\ 47rac J 

dependence, due to the finite optical depth of the envelope and the additional 
heating due to adiabatic contraction. 

4 APPROXIMATE METHOD OF COMPU¬ 
TATION OF THERMAL EVOLUTION 

The mathematical model describing the thermal structure of the protostellar cloud 
is fairly complex. The direct realization of the method used to numerically solve the 
corresponding system of equations in the multi-dimensional case involves the solu¬ 
tion of a number of additional problems of a technical nature. In particular, a highly 
sparse M x M matrix must be inverted, where M is the total number of cells in the 
computational grid. The classical methods for solving such problems (the method 
of Gauss, LU decomposition, etc.) are very inefficient. Complex iterative algo¬ 
rithms are used, such as multi-grid methods, the methods of Krylov, etc. [24]. The 
final choice of the method used in the multi-dimensional case should be determined 
by the desired compromise between accuracy, computational time, and complexity 
of the realization of the method. The use of adequate approximate methods in the 
multi-dimensional calculations is an attractive idea. We devised an approximate 
method for the solution of the system of equations describing the thermal structure 
of the protostellar cloud, which we call the “adaptive r-approximation”. In this 
method, the diffusion operator in the equation for the IR energy is replaced by the 
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local (coordinate dependent) algebraic operator 

We devised an approximate method for the solution of the system of equations 
describing the thermal structure of the protostellar cloud, which we call the “adap¬ 
tive r-approximation”. In this method, the diffusion operator in the equation for 
the IR energy is replaced by the local (coordinate dependent) algebraic operator 


V- 



Er -®bg 
T 


( 20 ) 


Here, E \> g is the energy density of the background IR radiation and r the charac¬ 
teristic diffusion time of the IR radiation. In essence, the adaptive r-approximation 
describes the release of IR radiation energy at a given point over a characteristic lo¬ 
cal timescale r, which is determined by the distribution of the radiation sources and 
the diffusion coefficient in the cloud. Generally speaking, various procedures can 
be used to calculate r in the algebraic operator. We chose r so that the stationary 
solution obtained using the r-approximation coincides with the exact solution. 

The coefficient r is computed as follows. Consider a stationary distribution 
of the gas and dust temperatures and the IR and UV energy. In this case, the 
subsystem of equations describing the thermal structure of the cloud can be reduced 
to the equation 


where 


V- 

Q 


—ve t ) 

= Q, 

(21) 

Or ) 


: Ag - r g 

- s. 

(22) 


Here, A g is the gas cooling function via the emission of molecular line radia¬ 
tion, T g the gas heating function via cosmic rays and photoelectrons, and S the 
dust heating function via external UV radiation. A bar denotes the steady-state 
distribution of the energy density. Solving this equation taking into account of the 
necessary boundary conditions yields an approximate distribution for the energy 
density of the IR radiation. The coefficient t can then be found using (20) and the 
corresponding approximate equation 


- E ' Ehs = Q. (23) 

T 

Technically, the adaptive r-approximation method can be implemented substan¬ 
tially more easily than the “full” method described above. Replacing the diffusion 
operator with an algebraic operator leads to a local nonlinear system of equations 
for determining T” +1 , T^ +1 , and T r n+1 . In turn, the finite-difference equation (21) 
is linear (<tr and Q depend on the previous time step) and its solution is appreciably 
simpler than the solution of the general system of equations. 

Consider the spherically symmetric case, assuming that the energy density of the 
IR radiation is equal to the background energy density at the outer boundary 
of the cloud, at r = R. At the inner boundary of the cloud, r = rq, the flux of 
IR radiation from the central protostar is specified to be F = T ce n. Integration of 
equations (21) and (23) yields 


r = 


1 

Q 



drirj 2 Q(ri) + rfF ceU 


Lri 


(24) 
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Note that the IR luminosity of the protostar is then L ce \\ = 47r?’i 2 i ? ce ii. In the 
case when cjr = const and Q = const, we have 


rjF ceU 


1 

R 


+ r R 


R 2 — r 2 


- n 


l 

R 


(25) 


For regions sufficiently far from the cloud boundary, we can assume r\ <C R and 
r < R. Therefore, we can obtain from (25) the approximate expression 


r 


^l^cellO-R . 1 d2 

-^ ■ 

Q r 6 


(26) 


Thus, the possibility of taking into account the additional source of IR ra¬ 
diation due to the central protostar is an important advantage of the adaptive 
r-approximation method. This makes it possible to use this method to calculate 
the thermal structure of a protostellar cloud in the multi-dimensional case with the 
constructed model for the protostar. 

The adaptive r-approximation method gives the exact solution (or close to it) 
in two important special cases. In the initial stages of the contraction of the pro¬ 
tostellar cloud, the medium is optically thin to IR radiation. This means that the 
timescale r is much shorter than the dynamical timescale, which can be estimated 
as the cloud free-fall timescale tff. As a result, the diffusion operator dominates, 
and we obtain nearly the steady-state energy density of the IR radiation, close to 
the background value Ey ig . Later, r increases proportional to the density in the 
central regions of the cloud, while dynamical timescale decreases with increasing 
density proportional to p Therefore, in late stages of the contraction, the sit¬ 
uation will become exactly the opposite: r will be much longer than the dynamical 
timescale. In this case, the diffusion operator is not important, and the energy 
density of the IR radiation will vary mainly due to the source terms. Therefore, we 
expect that the r-approximation will yield a solution close to the solution obtained 
using the rigorous diffusion operator. 

Figure 3 shows the distribution of the dust temperature Tj for various epochs 
in the evolution of a protostellar cloud obtained using the r-approximation (tau), 
Eddington approximation (edd), and diffusion approximation with a flux limiter 
(fid) (for a description of the Eddington method and the diffusion approximation 
with a flux limiter, see Section 2). The r-approximation and Eddington approxi¬ 
mation agree fairly well with the “fid” method in the prestellar phase, up to the 
time of the formation of the hydrostatic core. Significant differences between the 
distributions are observed in the accretion phase. The decrease in the dust tem¬ 
perature at 100-1000 AU is stronger in the Eddington approximation than in the 
“FLD” approach. This is associated with the ill-posed nature of the Eddington 
approximation for the streaming regime. The temperature obtained is higher in 
the r-approximation than in the Eddington approximation, but lower than in the 
“FLD” approximation. Although the adaptive r-approximation yields a significant 
error, it reproduces the behavior of the temperature fairly well, and is suitable for 
simple (preliminary) calculations of the evolution of protostellar clouds, as well as 
calculations of the evolution of optically thick protostellar disks. 
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Figure 3: Distributions of the dust temperature T d for various epochs during the evolution 
of a protostellar cloud, obtained using the r-approximation (tau), Eddington approxima¬ 
tion (edd), and the diffusion approximation with a flux limiter (fid). 








5 COMPARISON WITH OBSERVATIONS OF 
PRESTELLAR AND PROTOSTELLAR CORES 


The main aim of the thermal model presented here is a direct comparison of the 
results of the dynamical calculations with observations. The Herschel Space Ob¬ 
servatory has obtained rich observational data on various protostellar objects in 
recent years. In particular, the spatial and thermal structure of a number of dense 
molecular clouds, including both prestellar and protostellar cores, were studied in 
[25—27]. The distributions of the surface number density of hydrogen IV£ bs and 
the dust temperature averaged over the line of sight T° hs for a number of prestellar 
and protostellar cores are presented in [27]. It is natural for us is to attempt to 
reproduce these distributions using our model. Two circumstances must be taken 
into account when modeling the distributions of Afo and T d based on the theoret¬ 
ical distributions «h and T d . First, the temperature is not an additive quantity, 
and the averaging along the line of sight must be carried out consistently with the 
method used in [27] to reconstruct the distributions. Second, A^ bs and T^ hs were 
determined from observations of distributions of the radiation intensity with finite 
angular resolution, HPBW = 36 ,/ . 

Let us consider the radiation intensity I u in the optically thin approximation: 

Ah 

I u = m.HKv j B u (T d )dN , (27) 

o 

where k v [cm 2 /g] is the opacity per gram of gas, mu the atomic mass of 
hydrogen, N the surface number density of hydrogen along the line of sight, Afo 
the total surface number density of hydrogen, and B u the Planck function. In 
essence, it was assumed in [27] that 


/° bs = muKvB u (Ttm hs . (28) 

In other words, distributions T^ bs and A^ bs are found for which the spatial dis¬ 
tribution (28) satisfactorily describes the observed maps at all available frequencies. 
Since J° bs is a convolution of the true observed radiation intensity with the telescope 
beam, the distribution. T° hs d is a result of averaging not only along the line of 
sight, but also along the angular coordinate. 

To construct the theoretical distribution T^ onv , W e computed the distribution of 
the radiation intensity I u (x,y) using Eq. (27) and convolved it with the Gaussian 
telescope beam W(x 1 y)\ 


rconv 

In this formula, W(x,y ) = 


(x, y) = = J W(x - x',y - y')I„(x ', y')dx'dy'. 

1 


nH 2 


exp 


x 2 + y 2 \ x __ HPBW 
——~— , where H = - ; -. and 

H 2 J 2 v / hP2 


HPBW is the half-power beam width in radians. We compared T^ hs with the 
temperature . T^ onv found from the expression 


B u (tr v ) = 


rconv 


rconv 

' / H 


mu^Nfi 

Averaged over the telescope beam, the surface number density is 


(29) 
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Figure 4: Distributions for the prestellar core: hydrogen number density (upper left); 
gas, dust, and radiation temperature (upper right); surface density (lower left); and 
dust temperature weighted over the line of sight (lower right). In the lower panels, 
the model distributions are shown by the solid curves and the observed distributions 
from [27] by the bold dashed curves. The model distributions were derived by convolving 
the ideal distributions (shown by the thin dashed curves) with a Gaussian beam with 
HPBW = 36". 


A^ onv = = j W(x -x',y- y')N n {x ', y')dx'dy'. 

Note that the method for averaging the temperature described above does not 
require knowledge of the opacity k u in (27) and (30), since it cancels out in the 
final formula. However, the method does depend on the wavelength used in the 
Planck function. We used a wavelength of 160 mkm, close to the maximum of the 
observed radiation flux. We compared the model with the observed distributions 
averaged over the sample of sources shown in Fig. 8 of [27] by blue and red lines. 
When calculating the distributions N^ nv and T conv , we assumed that the cloud 
was located at a distance of 140 pc, and that the telescope beam was Gaussian 
with HPBW = 36". 
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Figure 5: Same as in Fig. 4 for protostcllar core. 


Table 1: Parameters of the prestellar cloud model 


Parameter 

Value 

Central concentration H 2 

4 x 10 4 cnr 3 

Radius 

5 x 10 4 AU 

Mass 

6 M 0 

Temperature of background UV radiation 

10 4 K 

Dilution of background UV radiation 

1.5 x 10~ 15 

Evolutionary time span 

330 ka 


We suppose that the mean observed distributions of the density and temper¬ 
ature for the prestellar cores correspond to cores close to hydrostatic and thermal 
equilibrium. Thus, the main parameters of the model are the initial central density, 
the radius of the cloud, the parameters of the external radiation field, and the evo¬ 
lution time. We selected the model parameters for the prestellar cores so that the 
value lV£ onv = 2 X 10 22 cnr 2 in the direction toward the cloud center coincided with 
the observed value, and the cloud radius was equal to the maximum value in the 
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distribution presented in [27], namely, R = 5 x 10 4 AU. The selected parameters 
of the external radiation field were such that Tj :onv = 15 K at the cloud boundary. 
The evolution time for the model was 1.3 f t fj; during this time, the central density 
doubled compared to the initial state. Note that the model parameters were not 
rigorously varied, and therefore are not strictly optimal, as we did not aim to ac¬ 
curately reproduce the observational data. Our goal was to determine whether our 
model was consistent with the available observational data. 

Figure 4 shows the distributions of the density and temperature in the model for 
the prestellar core, together with a comparison between the derived distributions 
N^ nv and Tj :onv and the observed distributions. The model parameters are given 
in the table. 

The results were obtained in the diffusion approximation with a flux limiter. 
The model and observed temperature distributions agree both qualitatively and 
quantitatively. The temperature is lower in the direction toward the center than 
in the direction toward the outer region of the cloud, as a consequence of external 
heating. At the same time, the surface density distributions differ significantly. 
The theoretical distribution has a more extended central plateau and steepe density 
decline toward the edge. At a distance of 2 x 10 4 AU, the model density is a factor of 
three higher than the observed value. One possible reason for this discrepancy in the 
observed and theoretical surface-density distributions is that the assumption that 
the prestellar cloud was initially in a quasi-equilibrium state was too crude. Indeed, 
according to modern gravo-turbulent concept of star formation [28], prestellar cores 
form in the turbulent interstellar medium. As a result of collisions of turbulent 
flows, gravitationally bound fragments form, whose density distributions may differ 
significantly from equilibrium distributions. 

We suppose that the mean observed density and temperature distributions for 
the protostellar cores correspond to clouds that are a product of the evolution of 
the prestellar cores and are in a phase of developed accretion. Therefore, the main 
parameter of the model is the evolution time. We choose this time to be t = 
3.5 tff, when the observed central surface number density Nff nv = 3 x 10 22 cnh’ 
is reproduced. Figure 5 shows model distributions of the density and temperature 
at this time together with the corresponding distributions of Ah and T and the 
observed profiles. 

At this time, the age of the star in the model is 340 000 yrs, its mass is 2.9 M 0 , 
the accretion rate is 6 x 10~ 6 Mq/ yr, the photospheric luminosity of the star is 
23 Lq, and the accretion luminosity is 45 L 0 . Overall, the theoretical and ob¬ 
served temperature distributions agree well. Note also that it is necessary to take 
the convolution into account when interpreting the observations, since the model 
temperature distribution unweighted by the telescope beam (thin solid curve) dif¬ 
fers strongly from the convolved distribution. At the same time, small differences 
in the distributions of the surface number density are observed in the outer parts of 
the cloud. This is obviously related to the initial conditions, similar to the situation 
with the prestellar cores. 

6 CONCLUSION 

We have presented a modification of amodel developed earlier for the calculation 
of the thermal structure of a collapsing protostellar cloud. The modified model is 
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intended for the calculation of the prestellar evolutionary stage of the cloud, as well 
as modeling of the accreting envelope after the formation of the protostar. The 
main improvements of the model are the introduction of the sink-cell formalism, 
calculation of the luminosity using a model for the evolution of the young star, 
and replacement of the Eddington approximation with the diffusion approximation 
with a flux limiter in the treatment of the radiative transfer. 

We have applied the model to the evolution of a protostellar cloud, right up to 
the formation of the first hydrostatic core and the subsequent accretion phase. In 
the inner region of the accreting envelope (r < 10 3 —10 4 AU), the temperature is 
governed by the radiation from the young star. In this zone, the temperatures of 
the gas and dust are equal and are well described by the distribution T oc r -1 / 2 . 
The thermal structure of the outer layers of the envelope (r > 10 3 —10 4 AU) is 
determined by interstellar radiation; the gas and dust temperatures rise toward the 
edge and may differ significantly. 

We have used our model to explain the distributions of the surface density and 
temperature averaged over the line of sight derived fromobservations of low mass 
prestellar and protostellar cores of molecular clouds obtained with the Herschel 
Space Observatory. The model and observed temperature distributions are in good 
agreement, for both the prestellar and protostellar cores. At the same time, the 
distributions of the surface density display some differences, probably due to the 
incorrect assumption that the cloud was initially in a quasi-hydrostatic state. The 
model can be adapted for multi-dimensional calculations, and can also be used to 
calculate the chemical structure of protostellar objects. 


The authors acknowledge Ralf Launhardt for providing observational data and 
B.M. Shustov for useful discussions. 
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